{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from math import sqrt, exp, erf\n",
    "import numpy as np\n",
    "\n",
    "def Sub_InitialField_SubLoop(i0, j0, k0, XI, YJ, ZK, model_EC, dx, dy, dz, subloop_x, subloop_y, L_subloop, Alpha):  \n",
    "    pi = 3.141592653589793\n",
    "    Permeability_vac = 1.256637061435917*10**-6\n",
    "    dt = np.empty((2))\n",
    "    EX_subloop = np.empty((XI, YJ+1, ZK+1, 2))\n",
    "    EY_subloop = np.empty((XI+1, YJ, ZK+1, 2))\n",
    "    EZ_subloop = np.empty((XI+1, YJ+1, ZK, 2))\n",
    "    HX_subloop = np.empty((XI+1, YJ, ZK, 2))\n",
    "    HY_subloop = np.empty((XI, YJ+1, ZK, 2))\n",
    "    HZ_subloop = np.empty((XI, YJ, ZK+1, 2))\n",
    "\n",
    "    initialfield_EC = model_EC[i0-1, j0-1, k0]\n",
    "    Txmm = L_subloop*L_subloop\n",
    "    t0_E = 1.13*Permeability_vac*model_EC[i0-1, j0-1, k0-1]*dz[k0-1]**2\n",
    "    dt[0] = Alpha*dz[k0-1] *         sqrt(Permeability_vac*model_EC[i0-1, j0-1, k0-1]*t0_E/6)\n",
    "    t0_H = t0_E+dt[0]/2\n",
    "\n",
    "    for k in range(1, ZK):\n",
    "        for j in range(1, YJ):\n",
    "            for i in range(0, XI):\n",
    "                coord_x = coor_x_calc(i, 1, i0, dx)+subloop_x\n",
    "                coord_y = coor_y_calc(j, 1, j0, dy)+subloop_y\n",
    "                coord_z = coor_z_calc(k, 1, k0, dz)\n",
    "                coord_r = sqrt(coord_x**2+coord_y**2+coord_z**2)\n",
    "                u = coord_r*sqrt(2*pi*initialfield_EC/(10**7*t0_E))\n",
    "                if u > 25:\n",
    "                    expu2 = float('inf')\n",
    "                else:\n",
    "                    expu2 = exp(u**2/2)\n",
    "                EX_subloop[i, j, k, 0] = sqrt(2/pi)*Txmm*coord_y*u**5/(4*pi*coord_r**5*initialfield_EC*expu2)\n",
    "    # Calculate subloop initial field Ey:\n",
    "    for k in range(1, ZK):\n",
    "        for j in range(0, YJ):\n",
    "            for i in range(1, XI):\n",
    "                coord_x = coor_x_calc(i, 2, i0, dx)+subloop_x\n",
    "                coord_y = coor_y_calc(j, 2, j0, dy)+subloop_y\n",
    "                coord_z = coor_z_calc(k, 2, k0, dz)\n",
    "                coord_r = sqrt(coord_x**2+coord_y**2+coord_z**2)\n",
    "                u = coord_r*sqrt(2*pi*initialfield_EC/(10**7*t0_E))\n",
    "                if u > 25:\n",
    "                    expu2 = float('inf')\n",
    "                else:\n",
    "                    expu2 = exp(u**2/2)\n",
    "                EY_subloop[i, j, k, 0] = sqrt(2/pi)*Txmm*coord_x*u**5/(4*pi*coord_r**5*initialfield_EC*expu2)\n",
    "\n",
    "    # Calculate subloop initial field Ez:\n",
    "    EZ_subloop[1:XI, 1:YJ, 0:ZK, 0] = 0\n",
    "    # Calculate subloop initial field Hx:\n",
    "    for k in range(0, ZK):\n",
    "        for j in range(0, YJ):\n",
    "            for i in range(0, XI+1):\n",
    "                coord_x = coor_x_calc(i, 4, i0, dx)+subloop_x\n",
    "                coord_y = coor_y_calc(j, 4, j0, dy)+subloop_y\n",
    "                coord_z = coor_z_calc(k, 4, k0, dz)\n",
    "                coord_r = sqrt(coord_x**2+coord_y**2+coord_z**2)\n",
    "                u = coord_r*sqrt(2*pi*initialfield_EC/(10**7*t0_H))\n",
    "                if u > 25:\n",
    "                    expu2 = float('inf')\n",
    "                else:\n",
    "                    expu2 = exp(u**2/2)\n",
    "                HX_subloop[i, j, k, 0] = Txmm*coord_z*coord_x/(4*pi*coord_r**5)*(3*erf(u)-sqrt(2/pi)*u*(3+u**2)/expu2)\n",
    "\n",
    "    # Calculate subloop initial field Hy:\n",
    "    for k in range(0, ZK):\n",
    "        for j in range(0, YJ+1):\n",
    "            for i in range(0, XI):\n",
    "                coord_x = coor_x_calc(i, 5, i0, dx)+subloop_x\n",
    "                coord_y = coor_y_calc(j, 5, j0, dy)+subloop_y\n",
    "                coord_z = coor_z_calc(k, 5, k0, dz)\n",
    "                coord_r = sqrt(coord_x**2+coord_y**2+coord_z**2)\n",
    "                u = coord_r*sqrt(2*pi*initialfield_EC/(10**7*t0_H))\n",
    "                if u > 25:\n",
    "                    expu2 = float('inf')\n",
    "                else:\n",
    "                    expu2 = exp(u**2/2)\n",
    "                HY_subloop[i, j, k, 0] = Txmm*coord_z*coord_y/(4*pi*coord_r**5)*(3*erf(u)-sqrt(2/pi)*u*(3+u**2)/expu2)\n",
    "\n",
    "    # Calculate subloop initial field Hz:\n",
    "    for k in range(0, ZK+1):\n",
    "        for j in range(0, YJ):\n",
    "            for i in range(0, XI):\n",
    "                coord_x = coor_x_calc(i, 6, i0, dx)+subloop_x\n",
    "                coord_y = coor_y_calc(j, 6, j0, dy)+subloop_y\n",
    "                coord_z = coor_z_calc(k, 6, k0, dz)\n",
    "                coord_r = sqrt(coord_x**2+coord_y**2+coord_z**2)\n",
    "                u = coord_r*sqrt(2*pi*initialfield_EC/(10**7*t0_H))\n",
    "                if u > 25:\n",
    "                    expu2 = float('inf')\n",
    "                else:\n",
    "                    expu2 = exp(u**2/2)\n",
    "                if coord_r < 0.001*dz[k0]:\n",
    "                    coord_r = 0.001*dz[k0] # Avoid singularity.\n",
    "                HZ_subloop[i, j, k, 0] = Txmm/(4*pi*coord_r**5)*((2*coord_z**2-sqrt(coord_x**2+coord_y**2))*erf(u)\n",
    "                       - (2*coord_z**2-sqrt(coord_x**2+coord_y**2)*(1+u**2))\n",
    "                       * sqrt(2/pi)/expu2)\n",
    "    return(EX_subloop, EY_subloop, EZ_subloop, HX_subloop, HY_subloop, HZ_subloop, t0_E, t0_H)\n",
    "\n",
    "def coor_x_calc(x, m,i0,dx):\n",
    "    coor_x=0\n",
    "    if (m == 2 or m == 3 or m == 4):\n",
    "        if (x <= i0-1):\n",
    "            for i in range(x, i0):\n",
    "                coor_x=coor_x+dx[i]\n",
    "            coor_x=coor_x-0.5*dx[i0-1]\n",
    "            coor_x=-coor_x\n",
    "        elif(x > i0-1):\n",
    "            for i in range(i0-1, x):\n",
    "                coor_x=coor_x+dx[i]\n",
    "            coor_x=coor_x-0.5*dx[i0-1]\n",
    "    else:\n",
    "        if (x <= i0-1):\n",
    "            for i in range(x, i0):\n",
    "                coor_x=coor_x+dx[i]\n",
    "            coor_x=coor_x-0.5*(dx[x]+dx[i0-1])\n",
    "            coor_x=-coor_x\n",
    "        elif(x > i0-1):\n",
    "            for i in range(i0-1, x+1):\n",
    "                coor_x=coor_x+dx[i]\n",
    "            coor_x=coor_x-0.5*(dx[i0-1]+dx[x])\n",
    "    return coor_x\n",
    "\n",
    "def coor_y_calc(y, m,j0,dy):\n",
    "    coor_y = 0\n",
    "    if (m == 1 or m == 3 or m == 5):\n",
    "        if (y <= j0-1):\n",
    "            for i in range(y,j0):\n",
    "                coor_y = coor_y+dy[i]\n",
    "            coor_y = coor_y-0.5*dy[j0-1]\n",
    "            coor_y = -coor_y\n",
    "        elif(y > j0-1):\n",
    "            for i in range(j0-1,y):\n",
    "                coor_y = coor_y+dy[i]\n",
    "            coor_y = coor_y-0.5*dy[j0-1]\n",
    "    else:\n",
    "        if (y <= j0-1):\n",
    "            for i in range(y,j0):\n",
    "                coor_y = coor_y+dy[i]\n",
    "            coor_y = coor_y-0.5*(dy[y]+dy[j0-1])\n",
    "            coor_y = -coor_y\n",
    "        elif(y > j0-1):\n",
    "            for i in range(j0-1,y+1):\n",
    "                coor_y = coor_y+dy[i]\n",
    "            coor_y = coor_y-0.5*(dy[j0-1]+dy[y])\n",
    "    return coor_y        \n",
    "\n",
    "def coor_z_calc(z, m,k0,dz):\n",
    "    coor_z = 0\n",
    "    if (m == 3 or m == 4 or m == 5):\n",
    "        if (z <= k0-1):\n",
    "            for i in range(z,k0):\n",
    "                coor_z = coor_z+dz[i]\n",
    "            coor_z = coor_z-0.5*dz[z]\n",
    "            coor_z = -coor_z\n",
    "        elif(z > k0-1):\n",
    "            for i in range(k0,z+1):\n",
    "                coor_z = coor_z+dz[i]\n",
    "            coor_z = coor_z-0.5*dz[z]\n",
    "    else:\n",
    "        if (z <= k0-1):\n",
    "            for i in range(z,k0):\n",
    "                coor_z = coor_z+dz[i]\n",
    "            coor_z = -coor_z\n",
    "        elif(z > k0-1):\n",
    "            for i in range(k0,z):\n",
    "                coor_z = coor_z+dz[i]      \n",
    "    return coor_z      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
